********************************************************************************
********************************************************************************
************ Results from home-specific cost-benefit analyses

********************************************************************************
******************** graphs incorporating social cost of carbon
clear 
use "$maindir/Results/Model_Outputs/cba_scc.dta"

**** identifying the percentiles of net present benefits
sort year20_npv
cumul year20_npv, gen(pctile_year20)
replace pctile_year20 = round(pctile_year20, 0.01)
tostring pctile_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv, gen(pctile_year30)
replace pctile_year30 = round(pctile_year30, 0.01)
tostring pctile_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv, gen(pctile_year10)
replace pctile_year10 = round(pctile_year10, 0.01)
tostring pctile_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv, gen(pctile_year40)
replace pctile_year40 = round(pctile_year40, 0.01)
tostring pctile_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv, gen(pctile_disc6)
replace pctile_disc6 = round(pctile_disc6, 0.01)
tostring pctile_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv, gen(pctile_disc0)
replace pctile_disc0 = round(pctile_disc0, 0.01)
tostring pctile_disc0, replace force format(%9.2f)

foreach y in year20 year30 year10 year40 disc6 disc0 {
sort `y'_npv
gen obsid_`y' = _n
labmask obsid_`y', values(pctile_`y')
}

gen year30_npv_1000 = year30_npv/1000
gen min95_year30_npv_1000 = min95_year30_npv/1000
gen max95_year30_npv_1000 = max95_year30_npv/1000


** labeling the percentiles
gsort -year30_npv
sum obsid_year30 if pctile_year30=="0.25" & pctile_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsid_year30 if pctile_year30=="0.10" & pctile_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsid_year30 if pctile_year30=="0.05" & pctile_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsid_year30 if pctile_year30=="0.01" & pctile_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsid_year30 if pctile_year30=="0.75" & pctile_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsid_year30 if pctile_year30=="0.90" & pctile_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsid_year30 if pctile_year30=="0.95" & pctile_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsid_year30 if pctile_year30=="0.99" & pctile_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsid_year30 if year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsid_year30 if pctile_year30=="0.50" & pctile_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med

sum year30_npv, detail
sca meannpv = round(`r(mean)', 0.01)
dis meannpv

sum year30_npv_1000, detail
sca meannpv1000 = `r(mean)'

count
sca numobs = r(N)

*** finally plotting the rank of NPV per home 
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsid_year30, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)8, labsize(small)) yscale(range(-12 9)) ///
	xline(`=scalar(obs_top25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-8 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75)) legend(off) note("Note: 95% confidence intervals based on bootstrapped standard errors.") ///
	yline(`=scalar(meannpv1000)', lcolor(red) lpattern(dash)) ///
	text(`=scalar(meannpv1000)+0.1' 2700 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(red) size(small) just(left) orient(vertical)) ///
	addplot((pcarrowi -8 `=scalar(obs_top25+50)' -8 `=scalar(obs_top25)+550', color(green)) || ///
	(pcarrowi 8 `=scalar(obs_bot25-50)' 8 `=scalar(obs_bot25)-550', color(green)))	
graph export "$maindir/Results/Figures/npb_homerank_fullsamp_scc.png", replace width(2500)



**** graphing the net present benefits - only homes with at least 12 months of data
sort year20_npv
cumul year20_npv if nmonths==12, gen(pctile12_year20)
replace pctile12_year20 = round(pctile12_year20, 0.01)
tostring pctile12_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv if nmonths==12, gen(pctile12_year30)
replace pctile12_year30 = round(pctile12_year30, 0.01)
tostring pctile12_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv if nmonths==12, gen(pctile12_year10)
replace pctile12_year10 = round(pctile12_year10, 0.01)
tostring pctile12_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv if nmonths==12, gen(pctile12_year40)
replace pctile12_year40 = round(pctile12_year40, 0.01)
tostring pctile12_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv if nmonths==12, gen(pctile12_disc6)
replace pctile12_disc6 = round(pctile12_disc6, 0.01)
tostring pctile12_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv if nmonths==12, gen(pctile12_disc0)
replace pctile12_disc0 = round(pctile12_disc0, 0.01)
tostring pctile12_disc0, replace force format(%9.2f)

gen aux = 1
foreach y in year20 year30 year10 year40 disc6 disc0 {
sort nmonths `y'_npv
gen obsid12_`y' = sum(aux) if nmonths==12
labmask obsid12_`y', values(pctile12_`y')
}


gsort nmonths -year30_npv
sum obsid12_year30 if pctile12_year30=="0.25" & pctile12_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsid12_year30 if pctile12_year30=="0.10" & pctile12_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsid12_year30 if pctile12_year30=="0.05" & pctile12_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsid12_year30 if pctile12_year30=="0.01" & pctile12_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsid12_year30 if pctile12_year30=="0.75" & pctile12_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsid12_year30 if pctile12_year30=="0.90" & pctile12_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsid12_year30 if pctile12_year30=="0.95" & pctile12_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsid12_year30 if pctile12_year30=="0.99" & pctile12_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsid12_year30 if nmonths==12 & year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsid12_year30 if pctile12_year30=="0.50" & pctile12_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med


* average baseline NPV
sum year30_npv if nmonths==12, detail
sca meannpv = round(`r(mean)', 0.01)
matrix ml_ate = (`r(mean)')
sca nobs = `r(N)'
dis meannpv

sum year30_npv_1000 if nmonths==12, detail
sca meannpv1000 = `r(mean)'

count if nmonths==12
sca numobs = r(N)

* calculate benefit-cost ratio
sum obsid12_year30 if obsid12_year30!=. & obsid12_year30[_n-1]==.
sca topobs = `r(mean)'
sum obsid12_year30 if _n==_N
sca botobs = `r(mean)'

gen toptemp1 = sum(year30_benefits) if pctile12_year30>="0.75" & pctile12_year30!="."
gen toptemp2 = sum(Real_nonHS_cost) if pctile12_year30>="0.75" & pctile12_year30!="."
gen bottemp1 = sum(year30_benefits) if pctile12_year30<="0.25" & pctile12_year30!="."
gen bottemp2 = sum(Real_nonHS_cost) if pctile12_year30<="0.25" & pctile12_year30!="."
gen cetemp1 = sum(year30_benefits) if obsid12_year30>=obs_zeronet & obsid12_year30!=.
gen cetemp2 = sum(Real_nonHS_cost) if obsid12_year30>=obs_zeronet & obsid12_year30!=.

sum toptemp1
sca ben_top25 = r(max)
sum toptemp2
sca cost_top25 = r(max)
sca bcr_top25 = ben_top25/cost_top25

sum bottemp1
sca ben_bot25 = r(max)
sum bottemp2
sca cost_bot25 = r(max)
sca bcr_bot25 = ben_bot25/cost_bot25

sum cetemp1
sca ben_ce = r(max)
sum cetemp2
sca cost_ce = r(max)
sca bcr_ce = ben_ce/cost_ce

di bcr_top25
di bcr_bot25
di bcr_ce

drop toptemp1 toptemp2 bottemp1 bottemp2 cetemp1 cetemp2

*** rank benefits per home - at least 12 months of data
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsid12_year30 if nmonths==12, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)12, labsize(small) glcolor(gray%10)) yscale(range(-12 12)) ///
	xline(`=scalar(obs_top25)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75)) ///
	yline(`=scalar(meannpv1000)', lcolor(gs3) lpattern(dash)) ///
	text(`=scalar(meannpv1000)-4' 2000 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(gs3) size(small) just(left) angle(45)) ///
	addplot((pcarrowi -8 `=scalar(obs_top25+25)' -8 `=scalar(obs_top25)+225', color(gs8)) ///
	(pcarrowi 8 `=scalar(obs_bot25-25)' 8 `=scalar(obs_bot25)-225', color(gs8)) ///
	(pcarrowi 9 `=scalar(obs_top25)' 9 `=scalar(topobs)', color(black)) ///
	(pcarrowi 9 `=scalar(topobs)' 9 `=scalar(obs_top25)', color(black)) ///
	(pcarrowi 9 `=scalar(botobs)' 9 `=scalar(obs_bot25)', color(black)) ///
	(pcarrowi 9 `=scalar(obs_bot25)' 9 `=scalar(botobs)', color(black)) ///
	(pcarrowi 11 `=scalar(obs_zeronet)' 11 `=scalar(topobs)', color(black)) ///
	(pcarrowi 11 `=scalar(topobs)' 11 `=scalar(obs_zeronet)', color(black)) ) ///
	text(9.1 `=scalar(obs_top10)-80' "Benefit-cost ratio = `=round(scalar(bcr_top25), 0.01)'", placement(n) color(black) size(vsmall) just(center)) ///
	text(9.1 `=scalar(obs_bot10)+120' "Benefit-cost ratio = 0`=round(scalar(bcr_bot25), 0.01)'", placement(n) color(black) size(vsmall) just(center)) ///
	text(11.1 `=scalar(obs_top10)-500' "Benefit-cost ratio = `=round(scalar(bcr_ce), 0.01)'", placement(n) color(black) size(vsmall) just(center))
graph export "$maindir/Results/Figures/npb_homerank_12month_scc.pdf", replace



***************************************** costs per ton of CO2 abated
* 0.000707  metric tons CO2/kWh of electricity
* 0.0549 metric tons CO2/Mcf of gas

* mmbtu savings to Mcf:
gen realized_savings_mcf = realized_savings_mmbtu*0.9756
* mmbtu savings to kWh:
gen realized_savings_kwh = realized_savings_mmbtu*293.07107

* converting energy savings to tons of CO2 equivalent
sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */
gen realized_savings_co2 = realized_savings_mcf*0.0549*(1-elec_contrib) + realized_savings_kwh*0.000707*elec_contrib


** cost-weighted average of lifetime of retrofits for each home
* cost share for each upgrade category (excluding H&S)
merge 1:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(tot_act* nonHS_cost Real_nonHS_cost) keep(3) nogen
drop if Real_nonHS_cost==.

foreach x in Attic AirCon AirSeal Baseload Door WtHtr Foundation General Window Furnace WallIns {
gen `x'_costshare = tot_act`x'/nonHS_cost
}
egen sharecheck = rowtotal(*costshare)

* based on WeatherWorks documentation: "Table SIR1 - Retrofit Service Lives"
gen avg_lifetime20 = 20*AirSeal_costshare + 25*WallIns_costshare + 25*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 25*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime20 = round(avg_lifetime20)
gen avg_lifetime_months20 = avg_lifetime20*12
replace avg_lifetime_months20 = round(avg_lifetime_months20)

* assuming all measures have a 10-year lifespan
gen avg_lifetime10 = 10*AirSeal_costshare + 10*WallIns_costshare + 10*Attic_costshare + ///
	10*Window_costshare + 10*Door_costshare + 10*Baseload_costshare + 10*AirCon_costshare + ///
	10*Furnace_costshare + 10*Foundation_costshare + 10*General_costshare + 10*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime10 = round(avg_lifetime10)
gen avg_lifetime_months10 = avg_lifetime10*12
replace avg_lifetime_months10 = round(avg_lifetime_months10)

* insulation measures with 50 years lifespan
gen avg_lifetime30 = 20*AirSeal_costshare + 50*WallIns_costshare + 50*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 50*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime30 = round(avg_lifetime30)
gen avg_lifetime_months30 = avg_lifetime30*12
replace avg_lifetime_months30 = round(avg_lifetime_months30)

* insulation measures with 80 years lifespan
gen avg_lifetime40 = 20*AirSeal_costshare + 80*WallIns_costshare + 80*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 80*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime40 = round(avg_lifetime40)
gen avg_lifetime_months40 = avg_lifetime40*12
replace avg_lifetime_months40 = round(avg_lifetime_months40)


***** co2 savings over 10 years
gsort nmonths -year10_npv
matrix costperco2 = (. ,. ,.)
gen tot_co2_abate = 0
forval i = 1/10 {
gen co2_abate`i' = realized_savings_co2
replace co2_abate`i' = 0 if round_lifetime10<`i' & co2_abate`i'!=.
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year10>="0.75" & pctile12_year10!="."
gen toptemp2 = sum(year10_npv) if pctile12_year10>="0.75" & pctile12_year10!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year10<="0.25" & pctile12_year10!="."
gen bottemp2 = sum(year10_npv) if pctile12_year10<="0.25" & pctile12_year10!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year10!=.
gen fulltemp2 = sum(year10_npv) if obsid12_year10!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 20 years
gsort nmonths -year20_npv
gen tot_co2_abate = 0
forval i = 1/25 {
gen co2_abate`i' = realized_savings_co2
replace co2_abate`i' = 0 if round_lifetime20<`i' & co2_abate`i'!=.
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year20>="0.75" & pctile12_year20!="."
gen toptemp2 = sum(year20_npv) if pctile12_year20>="0.75" & pctile12_year20!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year20<="0.25" & pctile12_year20!="."
gen bottemp2 = sum(year20_npv) if pctile12_year20<="0.25" & pctile12_year20!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year20!=.
gen fulltemp2 = sum(year20_npv) if obsid12_year20!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 30 years
gsort nmonths -year30_npv
gen tot_co2_abate = 0
forval i = 1/30 {
gen co2_abate`i' = realized_savings_co2
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year30>="0.75" & pctile12_year30!="."
gen toptemp2 = sum(year30_npv) if pctile12_year30>="0.75" & pctile12_year30!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year30<="0.25" & pctile12_year30!="."
gen bottemp2 = sum(year30_npv) if pctile12_year30<="0.25" & pctile12_year30!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year30!=.
gen fulltemp2 = sum(year30_npv) if obsid12_year30!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 40 years
gsort nmonths -year40_npv
gen tot_co2_abate = 0
forval i = 1/40 {
gen co2_abate`i' = realized_savings_co2
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year40>="0.75" & pctile12_year40!="."
gen toptemp2 = sum(year40_npv) if pctile12_year40>="0.75" & pctile12_year40!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year40<="0.25" & pctile12_year40!="."
gen bottemp2 = sum(year40_npv) if pctile12_year40<="0.25" & pctile12_year40!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year40!=.
gen fulltemp2 = sum(year40_npv) if obsid12_year40!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


matrix costperco2 = costperco2[2..., 1...]
matrix colnames costperco2 = "All Homes" "Top 25%" "Bottom 25%"
matrix rownames costperco2 = "10-year lifespans" "20-year lifespans" "30-year lifespans" "40-year lifespans"

matrix list costperco2 
esttab matrix(costperco2, fmt("%9.2f %9.2f %9.2f %9.2f")) ///
		using "$maindir/Results/Tables/costperco2_scc.tex", replace



**** graphing the net present benefits - prism sample restrictions
sort year30_npv
cumul year30_npv if prismrest==1, gen(pctileprism_year30)
replace pctileprism_year30 = round(pctileprism_year30, 0.01)
tostring pctileprism_year30, replace force format(%9.2f)

sort prismrest year30_npv
gen obsidprism = sum(aux) if prismrest==1
labmask obsidprism, values(pctileprism_year30)

gsort prismrest -year30_npv
sum obsidprism if pctileprism_year30=="0.25" & pctileprism_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsidprism if pctileprism_year30=="0.10" & pctileprism_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsidprism if pctileprism_year30=="0.05" & pctileprism_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsidprism if pctileprism_year30=="0.01" & pctileprism_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsidprism if pctileprism_year30=="0.75" & pctileprism_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsidprism if pctileprism_year30=="0.90" & pctileprism_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsidprism if pctileprism_year30=="0.95" & pctileprism_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsidprism if pctileprism_year30=="0.99" & pctileprism_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsidprism if prismrest==1 & year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsidprism if pctileprism_year30=="0.50" & pctileprism_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med

sum year30_npv if prismrest==1, detail
sca meannpv = round(`r(mean)', 0.01)
dis meannpv

sum year30_npv_1000 if prismrest==1, detail
sca meannpv1000 = `r(mean)'

count if prismrest==1
sca numobs = r(N)

*** rank benefits per home - prism
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsidprism if prismrest==1, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)8, labsize(small)) yscale(range(-12 9)) ///
	xline(`=scalar(obs_top25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-8 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75))  note("Note: 95% confidence intervals based on bootstrapped standard errors.") ///
	yline(`=scalar(meannpv1000)', lcolor(red) lpattern(dash)) ///
	text(`=scalar(meannpv1000)+0.15' 1700 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(red) size(small) just(left) orient(vertical)) ///
	addplot((pcarrowi -9 `=scalar(obs_top25+25)' -9 `=scalar(obs_top25)+225', color(green)) || ///
	(pcarrowi 9 `=scalar(obs_bot25-25)' 9 `=scalar(obs_bot25)-225', color(green)))
graph export "$maindir/Results/Figures/npb_homerank_prism_scc.png", replace width(2500)


***** table summarizing NPVs over the years
merge 1:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(ProgramYear) keep(3) nogen
keep if nmonths==12

********************************************************************************
********* net benefits by program year
gen post_arra = ProgramYear>2012

matrix summary_mat = (. , . , .)
forval i = 2009/2016 {
	qui: sum year30_npv if ProgramYear==`i'
	sca mean`i' = r(mean)
	sca sd`i' = r(sd)
	sca obs`i' = r(N)
	matrix tempmat = (mean`i', sd`i', obs`i')
	matrix summary_mat = (summary_mat \ tempmat)
}

qui: sum year30_npv if post_arra==0
sca meanpre = r(mean)
sca sdpre = r(sd)
sca obspre = r(N)

qui: sum year30_npv if post_arra==1
sca meanpost = r(mean)
sca sdpost = r(sd)
sca obspost = r(N)

matrix summary_mat = (summary_mat \ (., ., .) \ (meanpre, sdpre, obspre) \ (meanpost, sdpost, obspost))
matrix summary_mat = summary_mat[2..., 1...]

matrix colnames summary_mat = "Average" "Std. Dev." "Number of Homes" 
matrix rownames summary_mat = "PY 2009" "PY 2010" "PY 2011" "PY 2012" "PY 2013" "PY 2014" "PY 2015" "PY 2016" "" "PYs 2009-2012" "PYs 2013-2016"

esttab matrix(summary_mat, fmt(2 2 0)) ///
	using "$maindir/Results/Tables/cba_byPY_scc.tex", replace
		
		
		
		
		


********************************************************************************
********************** using retail energy prices
clear 
use "$maindir/Results/Model_Outputs/cba_retail.dta"

**** identifying the percentiles of net present benefits
sort year20_npv
cumul year20_npv, gen(pctile_year20)
replace pctile_year20 = round(pctile_year20, 0.01)
tostring pctile_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv, gen(pctile_year30)
replace pctile_year30 = round(pctile_year30, 0.01)
tostring pctile_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv, gen(pctile_year10)
replace pctile_year10 = round(pctile_year10, 0.01)
tostring pctile_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv, gen(pctile_year40)
replace pctile_year40 = round(pctile_year40, 0.01)
tostring pctile_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv, gen(pctile_disc6)
replace pctile_disc6 = round(pctile_disc6, 0.01)
tostring pctile_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv, gen(pctile_disc0)
replace pctile_disc0 = round(pctile_disc0, 0.01)
tostring pctile_disc0, replace force format(%9.2f)

foreach y in year20 year30 year10 year40 disc6 disc0 {
sort `y'_npv
gen obsid_`y' = _n
labmask obsid_`y', values(pctile_`y')
}

gen year30_npv_1000 = year30_npv/1000
gen min95_year30_npv_1000 = min95_year30_npv/1000
gen max95_year30_npv_1000 = max95_year30_npv/1000


** labeling the percentiles
gsort -year30_npv
sum obsid_year30 if pctile_year30=="0.25" & pctile_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsid_year30 if pctile_year30=="0.10" & pctile_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsid_year30 if pctile_year30=="0.05" & pctile_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsid_year30 if pctile_year30=="0.01" & pctile_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsid_year30 if pctile_year30=="0.75" & pctile_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsid_year30 if pctile_year30=="0.90" & pctile_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsid_year30 if pctile_year30=="0.95" & pctile_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsid_year30 if pctile_year30=="0.99" & pctile_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsid_year30 if year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsid_year30 if pctile_year30=="0.50" & pctile_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med

sum year30_npv, detail
sca meannpv = round(`r(mean)', 0.01)
dis meannpv

sum year30_npv_1000, detail
sca meannpv1000 = `r(mean)'

count
sca numobs = r(N)

*** finally plotting the rank of NPV per home 
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsid_year30, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)8, labsize(small)) yscale(range(-12 9)) ///
	xline(`=scalar(obs_top25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-8 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75)) legend(off) note("Note: 95% confidence intervals based on bootstrapped standard errors.") ///
	yline(`=scalar(meannpv1000)', lcolor(red) lpattern(dash)) ///
	text(`=scalar(meannpv1000)+0.1' 2700 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(red) size(small) just(left) orient(vertical)) ///
	addplot((pcarrowi -8 `=scalar(obs_top25+50)' -8 `=scalar(obs_top25)+550', color(green)) || ///
	(pcarrowi 8 `=scalar(obs_bot25-50)' 8 `=scalar(obs_bot25)-550', color(green)))	
graph export "$maindir/Results/Figures/npb_homerank_fullsamp_retail.png", replace width(2500)



**** graphing the net present benefits - only homes with at least 12 months of data
sort year20_npv
cumul year20_npv if nmonths==12, gen(pctile12_year20)
replace pctile12_year20 = round(pctile12_year20, 0.01)
tostring pctile12_year20, replace force format(%9.2f)

sort year30_npv
cumul year30_npv if nmonths==12, gen(pctile12_year30)
replace pctile12_year30 = round(pctile12_year30, 0.01)
tostring pctile12_year30, replace force format(%9.2f)

sort year10_npv
cumul year10_npv if nmonths==12, gen(pctile12_year10)
replace pctile12_year10 = round(pctile12_year10, 0.01)
tostring pctile12_year10, replace force format(%9.2f)

sort year40_npv
cumul year40_npv if nmonths==12, gen(pctile12_year40)
replace pctile12_year40 = round(pctile12_year40, 0.01)
tostring pctile12_year40, replace force format(%9.2f)

sort disc6_npv
cumul disc6_npv if nmonths==12, gen(pctile12_disc6)
replace pctile12_disc6 = round(pctile12_disc6, 0.01)
tostring pctile12_disc6, replace force format(%9.2f)

sort disc0_npv
cumul disc0_npv if nmonths==12, gen(pctile12_disc0)
replace pctile12_disc0 = round(pctile12_disc0, 0.01)
tostring pctile12_disc0, replace force format(%9.2f)

gen aux = 1
foreach y in year20 year30 year10 year40 disc6 disc0 {
sort nmonths `y'_npv
gen obsid12_`y' = sum(aux) if nmonths==12
labmask obsid12_`y', values(pctile12_`y')
}


gsort nmonths -year30_npv
sum obsid12_year30 if pctile12_year30=="0.25" & pctile12_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsid12_year30 if pctile12_year30=="0.10" & pctile12_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsid12_year30 if pctile12_year30=="0.05" & pctile12_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsid12_year30 if pctile12_year30=="0.01" & pctile12_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsid12_year30 if pctile12_year30=="0.75" & pctile12_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsid12_year30 if pctile12_year30=="0.90" & pctile12_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsid12_year30 if pctile12_year30=="0.95" & pctile12_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsid12_year30 if pctile12_year30=="0.99" & pctile12_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsid12_year30 if nmonths==12 & year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsid12_year30 if pctile12_year30=="0.50" & pctile12_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med


* average baseline NPV
sum year30_npv if nmonths==12, detail
sca meannpv = round(`r(mean)', 0.01)
matrix ml_ate = (`r(mean)')
sca nobs = `r(N)'
dis meannpv

sum year30_npv_1000 if nmonths==12, detail
sca meannpv1000 = `r(mean)'

count if nmonths==12
sca numobs = r(N)

* calculate benefit-cost ratio
sum obsid12_year30 if obsid12_year30!=. & obsid12_year30[_n-1]==.
sca topobs = `r(mean)'
sum obsid12_year30 if _n==_N
sca botobs = `r(mean)'

gen toptemp1 = sum(year30_benefits) if pctile12_year30>="0.75" & pctile12_year30!="."
gen toptemp2 = sum(Real_nonHS_cost) if pctile12_year30>="0.75" & pctile12_year30!="."
gen bottemp1 = sum(year30_benefits) if pctile12_year30<="0.25" & pctile12_year30!="."
gen bottemp2 = sum(Real_nonHS_cost) if pctile12_year30<="0.25" & pctile12_year30!="."
gen cetemp1 = sum(year30_benefits) if obsid12_year30>=obs_zeronet & obsid12_year30!=.
gen cetemp2 = sum(Real_nonHS_cost) if obsid12_year30>=obs_zeronet & obsid12_year30!=.

sum toptemp1
sca ben_top25 = r(max)
sum toptemp2
sca cost_top25 = r(max)
sca bcr_top25 = ben_top25/cost_top25

sum bottemp1
sca ben_bot25 = r(max)
sum bottemp2
sca cost_bot25 = r(max)
sca bcr_bot25 = ben_bot25/cost_bot25

sum cetemp1
sca ben_ce = r(max)
sum cetemp2
sca cost_ce = r(max)
sca bcr_ce = ben_ce/cost_ce

di bcr_top25
di bcr_bot25
di bcr_ce

drop toptemp1 toptemp2 bottemp1 bottemp2 cetemp1 cetemp2

*** rank benefits per home - at least 12 months of data
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsid12_year30 if nmonths==12, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)12, labsize(small) glcolor(gray%10)) yscale(range(-12 12)) ///
	xline(`=scalar(obs_top25)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(gs8) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(black) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75)) ///
	yline(`=scalar(meannpv1000)', lcolor(gs3) lpattern(dash)) ///
	text(`=scalar(meannpv1000)+4' 2000 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(gs3) size(small) just(left) orient(horizontal)) ///
	addplot((pcarrowi -8 `=scalar(obs_top25+25)' -8 `=scalar(obs_top25)+225', color(gs8)) ///
	(pcarrowi 8 `=scalar(obs_bot25-25)' 8 `=scalar(obs_bot25)-225', color(gs8)) ///
	(pcarrowi 9 `=scalar(obs_top25)' 9 `=scalar(topobs)', color(black)) ///
	(pcarrowi 9 `=scalar(topobs)' 9 `=scalar(obs_top25)', color(black)) ///
	(pcarrowi 9 `=scalar(botobs)' 9 `=scalar(obs_bot25)', color(black)) ///
	(pcarrowi 9 `=scalar(obs_bot25)' 9 `=scalar(botobs)', color(black)) ///
	(pcarrowi 11 `=scalar(obs_zeronet)' 11 `=scalar(topobs)', color(black)) ///
	(pcarrowi 11 `=scalar(topobs)' 11 `=scalar(obs_zeronet)', color(black)) ) ///
	text(9.1 `=scalar(obs_top10)-80' "Benefit-cost ratio = `=round(scalar(bcr_top25), 0.01)'", placement(n) color(black) size(vsmall) just(center)) ///
	text(9.1 `=scalar(obs_bot10)+120' "Benefit-cost ratio = 0`=round(scalar(bcr_bot25), 0.01)'", placement(n) color(black) size(vsmall) just(center)) ///
	text(11.1 `=scalar(obs_top10)-800' "Benefit-cost ratio = `=round(scalar(bcr_ce), 0.01)'", placement(n) color(black) size(vsmall) just(center))
graph export "$maindir/Results/Figures/npb_homerank_12month_retail.pdf", replace



***************************************** costs per ton of CO2 abated
* 0.000707  metric tons CO2/kWh of electricity
* 0.0549 metric tons CO2/Mcf of gas

* mmbtu savings to Mcf:
gen realized_savings_mcf = realized_savings_mmbtu*0.9756
* mmbtu savings to kWh:
gen realized_savings_kwh = realized_savings_mmbtu*293.07107

* converting energy savings to tons of CO2 equivalent
sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */
gen realized_savings_co2 = realized_savings_mcf*0.0549*(1-elec_contrib) + realized_savings_kwh*0.000707*elec_contrib


** cost-weighted average of lifetime of retrofits for each home
* cost share for each upgrade category (excluding H&S)
merge 1:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(tot_act* nonHS_cost Real_nonHS_cost) keep(3) nogen
drop if Real_nonHS_cost==.

foreach x in Attic AirCon AirSeal Baseload Door WtHtr Foundation General Window Furnace WallIns {
gen `x'_costshare = tot_act`x'/nonHS_cost
}
egen sharecheck = rowtotal(*costshare)

* based on WeatherWorks documentation: "Table SIR1 - Retrofit Service Lives"
gen avg_lifetime20 = 20*AirSeal_costshare + 25*WallIns_costshare + 25*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 25*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime20 = round(avg_lifetime20)
gen avg_lifetime_months20 = avg_lifetime20*12
replace avg_lifetime_months20 = round(avg_lifetime_months20)

* assuming all measures have a 10-year lifespan
gen avg_lifetime10 = 10*AirSeal_costshare + 10*WallIns_costshare + 10*Attic_costshare + ///
	10*Window_costshare + 10*Door_costshare + 10*Baseload_costshare + 10*AirCon_costshare + ///
	10*Furnace_costshare + 10*Foundation_costshare + 10*General_costshare + 10*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime10 = round(avg_lifetime10)
gen avg_lifetime_months10 = avg_lifetime10*12
replace avg_lifetime_months10 = round(avg_lifetime_months10)

* insulation measures with 50 years lifespan
gen avg_lifetime30 = 20*AirSeal_costshare + 50*WallIns_costshare + 50*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 50*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime30 = round(avg_lifetime30)
gen avg_lifetime_months30 = avg_lifetime30*12
replace avg_lifetime_months30 = round(avg_lifetime_months30)

* insulation measures with 80 years lifespan
gen avg_lifetime40 = 20*AirSeal_costshare + 80*WallIns_costshare + 80*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 80*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime40 = round(avg_lifetime40)
gen avg_lifetime_months40 = avg_lifetime40*12
replace avg_lifetime_months40 = round(avg_lifetime_months40)


***** co2 savings over 10 years
gsort nmonths -year10_npv
matrix costperco2 = (. ,. ,.)
gen tot_co2_abate = 0
forval i = 1/10 {
gen co2_abate`i' = realized_savings_co2
replace co2_abate`i' = 0 if round_lifetime10<`i' & co2_abate`i'!=.
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year10>="0.75" & pctile12_year10!="."
gen toptemp2 = sum(year10_npv) if pctile12_year10>="0.75" & pctile12_year10!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year10<="0.25" & pctile12_year10!="."
gen bottemp2 = sum(year10_npv) if pctile12_year10<="0.25" & pctile12_year10!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year10!=.
gen fulltemp2 = sum(year10_npv) if obsid12_year10!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 20 years
gsort nmonths -year20_npv
gen tot_co2_abate = 0
forval i = 1/25 {
gen co2_abate`i' = realized_savings_co2
replace co2_abate`i' = 0 if round_lifetime20<`i' & co2_abate`i'!=.
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year20>="0.75" & pctile12_year20!="."
gen toptemp2 = sum(year20_npv) if pctile12_year20>="0.75" & pctile12_year20!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year20<="0.25" & pctile12_year20!="."
gen bottemp2 = sum(year20_npv) if pctile12_year20<="0.25" & pctile12_year20!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year20!=.
gen fulltemp2 = sum(year20_npv) if obsid12_year20!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 30 years
gsort nmonths -year30_npv
gen tot_co2_abate = 0
forval i = 1/30 {
gen co2_abate`i' = realized_savings_co2
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year30>="0.75" & pctile12_year30!="."
gen toptemp2 = sum(year30_npv) if pctile12_year30>="0.75" & pctile12_year30!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year30<="0.25" & pctile12_year30!="."
gen bottemp2 = sum(year30_npv) if pctile12_year30<="0.25" & pctile12_year30!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year30!=.
gen fulltemp2 = sum(year30_npv) if obsid12_year30!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


***** co2 savings over 40 years
gsort nmonths -year40_npv
gen tot_co2_abate = 0
forval i = 1/40 {
gen co2_abate`i' = realized_savings_co2
replace tot_co2_abate = tot_co2_abate + co2_abate`i'
drop co2_abate`i'
}

gen toptemp1 = sum(tot_co2_abate) if pctile12_year40>="0.75" & pctile12_year40!="."
gen toptemp2 = sum(year40_npv) if pctile12_year40>="0.75" & pctile12_year40!="."
gen bottemp1 = sum(tot_co2_abate) if pctile12_year40<="0.25" & pctile12_year40!="."
gen bottemp2 = sum(year40_npv) if pctile12_year40<="0.25" & pctile12_year40!="."
gen fulltemp1 = sum(tot_co2_abate) if obsid12_year40!=.
gen fulltemp2 = sum(year40_npv) if obsid12_year40!=.

sum toptemp1 if toptemp1[_n]!=. & toptemp1[_n+1]==.
sca ben_top25 = r(mean)
sum toptemp2  if toptemp2[_n]!=. & toptemp2[_n+1]==.
sca cost_top25 = r(mean)
sca co2_top25 = cost_top25/ben_top25

sum bottemp1 if _n==_N
sca ben_bot25 = r(mean)
sum bottemp2 if _n==_N
sca cost_bot25 = r(mean)
sca co2_bot25 = cost_bot25/ben_bot25

sum fulltemp1 if _n==_N
sca ben_ce = r(mean)
sum fulltemp2 if _n==_N
sca cost_ce = r(mean)
sca co2_ce = cost_ce/ben_ce

di co2_top25
di co2_bot25
di co2_ce

matrix costperco2 = (costperco2 \ co2_ce, co2_top25, co2_bot25)

drop toptemp1 toptemp2 bottemp1 bottemp2 fulltemp1 fulltemp2 tot_co2_abate


matrix costperco2 = costperco2[2..., 1...]
matrix colnames costperco2 = "All Homes" "Top 25%" "Bottom 25%"
matrix rownames costperco2 = "10-year lifespans" "20-year lifespans" "30-year lifespans" "40-year lifespans"

matrix list costperco2 
esttab matrix(costperco2, fmt("%9.2f %9.2f %9.2f %9.2f")) ///
		using "$maindir/Results/Tables/costperco2_retail.tex", replace



**** graphing the net present benefits - prism sample restrictions
sort year30_npv
cumul year30_npv if prismrest==1, gen(pctileprism_year30)
replace pctileprism_year30 = round(pctileprism_year30, 0.01)
tostring pctileprism_year30, replace force format(%9.2f)

sort prismrest year30_npv
gen obsidprism = sum(aux) if prismrest==1
labmask obsidprism, values(pctileprism_year30)

gsort prismrest -year30_npv
sum obsidprism if pctileprism_year30=="0.25" & pctileprism_year30[_n-1]=="0.26"
sca obs_bot25 = `r(mean)'
dis obs_bot25

sum obsidprism if pctileprism_year30=="0.10" & pctileprism_year30[_n-1]=="0.11"
sca obs_bot10 = `r(mean)'
dis obs_bot10

sum obsidprism if pctileprism_year30=="0.05" & pctileprism_year30[_n-1]=="0.06"
sca obs_bot5 = `r(mean)'
dis obs_bot5

sum obsidprism if pctileprism_year30=="0.01" & pctileprism_year30[_n-1]=="0.02"
sca obs_bot1 = `r(mean)'
dis obs_bot1

sum obsidprism if pctileprism_year30=="0.75" & pctileprism_year30[_n-1]=="0.76"
sca obs_top25 = `r(mean)'
dis obs_top25

sum obsidprism if pctileprism_year30=="0.90" & pctileprism_year30[_n-1]=="0.91"
sca obs_top10 = `r(mean)'
dis obs_top10

sum obsidprism if pctileprism_year30=="0.95" & pctileprism_year30[_n-1]=="0.96"
sca obs_top5 = `r(mean)'
dis obs_top5

sum obsidprism if pctileprism_year30=="0.99" & pctileprism_year30[_n-1]=="1.00"
sca obs_top1 = `r(mean)'
dis obs_top1

sum obsidprism if prismrest==1 & year30_npv<=0 & year30_npv[_n-1]>0
sca obs_zeronet = `r(mean)'
dis obs_zeronet

sum obsidprism if pctileprism_year30=="0.50" & pctileprism_year30[_n-1]=="0.51"
sca obs_med = `r(mean)'
dis obs_med

sum year30_npv if prismrest==1, detail
sca meannpv = round(`r(mean)', 0.01)
dis meannpv

sum year30_npv_1000 if prismrest==1, detail
sca meannpv1000 = `r(mean)'

count if prismrest==1
sca numobs = r(N)

*** rank benefits per home - prism
eclplot year30_npv_1000 min95_year30_npv_1000 max95_year30_npv_1000 obsidprism if prismrest==1, ///
	eplottype(line) rplottype(rarea) ///
	estopts(lcolor(black) lwidth(thick))  ///
	ciopts(fcolor(gray%75) lwidth(none))  ///
	xtitle("Home Percentile Rank: best homes {&rarr}") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("IHWAP Net Present Benefits ($1000)") ///
	xlabel(`=scalar(obs_bot25)', labsize(small) angle(45) valuelabel) ///
	xlabel(`=scalar(obs_med)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top25)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_bot10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_top10)', labsize(small) angle(45) valuelabel add) ///
	xlabel(`=scalar(obs_zeronet)', labsize(small) angle(45) valuelabel add) ///
	ylabel(-12(2)8, labsize(small)) yscale(range(-12 9)) ///
	xline(`=scalar(obs_top25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top25)+0.07' "Top 25% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top10)+0.07' "Top 10% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top5)+0.07' "Top 5% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_top1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(-4.2 `=scalar(obs_top1)+0.07' "Top 1% homes", placement(e) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot25)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot25)-0.07' "Bottom 25% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot10)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot10)-0.07' "Bottom 10% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot5)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot5)-0.07' "Bottom 5% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_bot1)', lcolor(green) lpattern(shortdash_dot)) ///
	text(4.2 `=scalar(obs_bot1)-0.07' "Bottom 1% homes", placement(w) color(green) size(small) just(left) orient(vertical)) ///
	xline(`=scalar(obs_zeronet)', lcolor(black) lpattern(shortdash_dot)) ///
	text(-8 `=scalar(obs_zeronet)-0.07' "MB = MC", placement(e) color(black) size(small) just(left) orient(vertical)) ///
	yline(0, lcolor(gray%75))  note("Note: 95% confidence intervals based on bootstrapped standard errors.") ///
	yline(`=scalar(meannpv1000)', lcolor(red) lpattern(dash)) ///
	text(`=scalar(meannpv1000)+0.15' 1700 "Avg. Benefits = US$ `=round(scalar(meannpv))'", placement(n) color(red) size(small) just(left) orient(vertical)) ///
	addplot((pcarrowi -9 `=scalar(obs_top25+25)' -9 `=scalar(obs_top25)+225', color(green)) || ///
	(pcarrowi 9 `=scalar(obs_bot25-25)' 9 `=scalar(obs_bot25)-225', color(green)))
graph export "$maindir/Results/Figures/npb_homerank_prism_retail.png", replace width(5000)



********************************************************************************
***** table summarizing NPVs over the years
merge 1:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(ProgramYear) keep(3) nogen
keep if nmonths==12

********************************************************************************
********* net benefits by program year
gen post_arra = ProgramYear>2012

matrix summary_mat = (. , . , .)
forval i = 2009/2016 {
	qui: sum year30_npv if ProgramYear==`i'
	sca mean`i' = r(mean)
	sca sd`i' = r(sd)
	sca obs`i' = r(N)
	matrix tempmat = (mean`i', sd`i', obs`i')
	matrix summary_mat = (summary_mat \ tempmat)
}

qui: sum year30_npv if post_arra==0
sca meanpre = r(mean)
sca sdpre = r(sd)
sca obspre = r(N)

qui: sum year30_npv if post_arra==1
sca meanpost = r(mean)
sca sdpost = r(sd)
sca obspost = r(N)

matrix summary_mat = (summary_mat \ (., ., .) \ (meanpre, sdpre, obspre) \ (meanpost, sdpost, obspost))
matrix summary_mat = summary_mat[2..., 1...]

matrix colnames summary_mat = "Average" "Std. Dev." "Number of Homes" 
matrix rownames summary_mat = "PY 2009" "PY 2010" "PY 2011" "PY 2012" "PY 2013" "PY 2014" "PY 2015" "PY 2016" "" "PYs 2009-2012" "PYs 2013-2016"

esttab matrix(summary_mat, fmt(2 2 0)) ///
	using "$maindir/Results/Tables/cba_byPY_retail.tex", replace
